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Abstract: We compute the fully differential decay rate of the Standard Model Higgs boson into b- 
quarks at next-to-next-to-leading order (NNLO) accuracy in a s . We employ a general subtraction 
scheme developed for the calculation of higher order perturbative corrections to QCD jet cross 
sections, which is based on the universal infrared factorization properties of QCD squared matrix 
elements. We show that the subtractions render the various contributions to the NNLO correction 
finite. In particular, we demonstrate analytically that the sum of integrated subtraction terms 
correctly reproduces the infrared poles of the two-loop double virtual contribution to this process. 
We present illustrative differential distributions obtained by implementing the method in a parton 
level Monte Carlo program. The basic ingredients of our subtraction scheme, used here for the first 
time to compute a physical observable, are universal and can be employed for the computation of 
more involved processes. 
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1 Introduction 

In run I, the ATLAS and CMS collaborations of the Large Hadron Collider (LHC) discovered a new 
particle [1, 2] with quantum numbers corresponding to those of the Higgs boson in the Standard 
Model (SM) within the experimental accuracy of the measurements [3-6]. Thus by now it is widely 
accepted that the new particle is the Higgs boson of the SM. Nevertheless, further more precise 
measurements are being prepared for the upcoming run II. In particular, a lot of emphasis is put 
on the precise determination of the couplings of the Higgs boson to the heavy fermions to check 
whether the fermion masses are consistent with fermion mass generation in the SM. 

Since the b-quark is quite light (its mass is only about 2 % of the vacuum expectation value of 
the Higgs field), the rate of associated production of a b-quark pair with a Higgs boson is rather 
low. This fact, together with the overwhelming number of background events coming from direct 
QCD b-quark pair production makes the determination of the b-quark Yukawa coupling through 
H bb production impossible. A better option that gives direct access to the Hbb Yukawa coupling 
is to measure the H —» bb decay in the associated production of a Higgs boson with a W or a Z 
boson in a boosted or semi-boosted regime [7]. In this scenario it is possible to use the kinematic 
and topological properties of the final states to isolate the H —> bb decay. In this respect, first 
measurements have been performed by the CMS [8] and ATLAS [9] collaborations. 
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Such search strategies may be aided by accurate modeling of QCD radiation in the H —> bb 
decay, which motivates the computation of the fully differential decay rate at next-to-next-to- 
leading order (NNLO) accuracy in QCD perturbation theory. Computing fully differential cross 
sections and decay rates at NNLO turns out to be rather involved, however the last decade has 
witnessed substantial development [10-41] leading to a number of differential results for specific 
processes [42-76]. 

The first computation of the fully differential decay rate of the SM Higgs boson into b-quarks 
at NNLO accuracy was published in ref. [47] . That computation was performed with the method of 
sector decomposition based on non-linear mappings [13]. Here we offer a different approach based 
on the numerical implementation of the general subtraction scheme developed in a series of papers 
for the computation of QCD jet cross sections at NNLO accuracy [31-41]. This method, which 
is used for the first time in this paper to compute a physical observable at NNLO, employs the 
universal infrared factorization of QCD squared matrix elements to define local subtraction terms 
for regulating the singularities emerging in unresolved real radiation. 

Specifically, we can write the NNLO correction to the cross section of a generic m-jet process 
as a sum of three contributions, the tree level double real radiation, the one-loop plus a single 
radiation, and the two-loop double virtual terms of the basic process under consideration, 


^nnlo _ 
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The subscripts on the integral signs are simply reminders that the integration is over the phase space 
of n = to, to + 1 or m + 2 final state particles. Above J„ denotes the value of some infrared-safe 
observable J evaluated on an n parton final state. 

The right-hand sides of eqs. (1.3) and (1.4) are integrable in four dimensions by construction 
[31-34], while the integrability of eq. (1.5) in four dimensions is ensured by the Kinoshita-Lee- 
Nauenberg (KLN) theorem on infrared-safe quantities, provided that our subtraction scheme is well 
defined. 

The counterterms which contribute to d++ 2 ° and t° dcr^^ 0 were introduced in refs. [33] 
and [34]. The integration of the real-virtual counterterms (the last two terms of eq. (1.5)) was 
performed in refs. [35, 36, 38]. The integral of the iterated single unresolved counterterm (the third 
term of eq. (1.5)) was computed in ref. [39]. The integration of the collinear-type contributions to the 
double unresolved counterterm (the second term of eq. (1.5)) was performed in ref. [40]. The soft- 
type contributions to the same counterterm were presented in ref. [41]. Most of these results were 
given as expansions in e whose coefficients were computed numerically. Here we present the relevant 
integrals with pole coefficients evaluated analytically, while the finite parts are given numerically. 
The final test on the consistency of our subtraction scheme is then to verify that eq. (1.5) is free of 
singularities, as prescribed by the KLN theorem. In this paper, we perform that check analytically 
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for the first time by computing the fully differential decay rate 1 of the Higgs boson into b-quarks 
at NNLO. 

The present work is the first physical application of this method, therefore in order to facilitate 
reading we present the full computation as implemented in a parton level Monte Carlo program in 
detail. As usual in such codes, the jet function J is computed from generated momenta in d = 4 
dimensions, therefore, the implementation of any infrared-safe physical quantity is straightforward 
as demonstrated here. 

The paper is organized as follows: in section 2, the notation and conventions we use are in¬ 
troduced; in sections 3 and 4, we show the decay width at leading order and next-to-leading order 
(NLO) accuracy in a s ; in section 5, we display the counterterms and the insertion operators which 
are necessary to define the double real (1.3) and the real-virtual (1.4) contributions to the decay 
width, and we show that the double virtual contribution (1.5) is free of singularities; in section 6, 
we show a selection of illustrative results; we draw our conclusions in section 7. The two appendices 
provide details on the matrix elements we use, as well as on the insertion operator used in the NLO 
computation. 


2 Notation 

We consider the partial decay width of the Higgs boson into a b-quark pair, for any 

infrared-safe observable J. Through NNLO in QCD, this decay width receives contributions from 
the following partonic subprocesses: 


LO 

H(ph) -t b(pi) + b(p 2 ) 

tree level 

NLO 

H(ph) -> b(pi) + b(p 2 ) + g{ps) 

tree level 


H(ph) -t b(pi) + b(p 2 ) 

one-loop 

NNLO 

H(ph) -t b(pi) + b (p 2 ) + g{p 3 ) + g(p4) 

tree level 


H(ph ) -t b(pi) + b(p 2 ) + q(P3) + q(pa) 

tree level 


H(ph ) -t b(pi) + b(p 2 ) +b(p 3 ) + b(p 4 ) 

tree level 


H(ph) -t b(pi) + b(p 2 ) + g{p 3 ) 

one-loop 


H(pn) -t b(pi) + b (p 2 ) 

two-loop 


where we show also the four-momenta of the particles in parentheses. We report the matrix elements 
corresponding to all subprocesses up to the required loop level in appendix A. 

We use the colour and spin space notation of ref. [77] where the matrix element for a given 
subprocess, |A4 n ), is a vector in color and spin space, normalized such that the squared matrix 
element summed over colours and spins is given by 

\MJ 2 = (M n \M n ), (2.1) 

where n is the number of particles in the final state. The matrix element has the following formal 
loop expansion 

\M n ) = \M^) + \M^) + \M^) + ... , 

with the dots denoting higher-loop contributions. We will always consider matrix elements com¬ 
puted in conventional dimensional regularization (CDR) with MS subtraction. We will also use the 

1 In eqs. (1.1)—(1.5) we presented the basic structure of our subtraction scheme for computing a generic cross 

section, however our method applies equally to decay rates, as spelled out in detail in sections 3—5. 
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following ® product notation to indicate the insertion of colour charge operators between 
and \M^): 


{M^ l} \M (e2) ) ® T t -T k = (M (tl) \T r T k |aM } ), 
{M^\M^)®{T i -T k ,T j -Ti} = {M { ^\{T i -T k ,T j -Ti}\M^). 


( 2 . 2 ) 


We use the customary normalization of Tr = 1/2 for the colour-charge operators, thus the quadratic 
Casimirs are C A = 2T R N C = N c in the adjoint and C F = Tr(N 2 — 1 )/(N c ) = (N 2 — 1)/(2 N c ) in 
the fundamental representation, where N c = 3 is the number of colours. 

The b-quark mass is much smaller than the scale of the problem that is the Higgs boson mass, 
therefore, we treat the b-quarks as massless, both in the matrix elements and phase space integrals, 
retaining the b-quark mass only in the Yukawa coupling. We neglect the t-quark throughout and 
consider n f = 5 light quark flavours. 

In QCD the renormalized amplitudes are obtained from the unrenormalized ones by replacing 
the bare couplings y^ and af with their renormalized counterparts evaluated at the renormalization 
scale /I 
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where 


HC'a 4n f T R 
3 3 


(2.5) 


and = (47r) e exp(—gyp) corresponds to MS subtraction. Although the factor (47r) £ exp(— c7e) 
is often abbreviated as S e in the literature, we reserve the latter to denote 


, = (4t r) £ 

£ r(i- e )' 


( 2 . 6 ) 


On the right-hand side, yb = 2/b(M) an d a s = (m) are th e dimensionless renormalized couplings 

in the MS scheme evaluated at the renormalization scale p. 

The n particle massless phase space measure reads 


d (j>n (Q ) — d<p n (pi, . . . , p n , Q ) 
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d d Pi 
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Throughout the paper, we will use yi k to denote twice the dot-product of two momenta, scaled by 
the total momentum squared Q 2 . For example, 


_ 2 pi ■ p k 

y ik - Q2 


and 


_ 2,pi ■ Q 
VlQ ~ ~Q 2 ~ 


We also introduce the combination 


Yik,Q = 


Vik 

ViQUkQ 


for later convenience. 


( 2 . 8 ) 


(2.9) 
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3 Leading order 

Let us denote the Born differential decay rate by, 


dr 2 = l-^b&l 2 • ( 3 ' 1 ) 

Then the leading order decay width is, 

r B [j] = dr B j 2 = j d Mm 2 H ) \m { 2\ 2 Ji ■ (3.2) 

Here J is an infrared-safe observable whose value evaluated on a kinematic configuration with two 
partons is J 2 . For the inclusive decay width (J = 1) at leading order we have 

T LO = T b [J = 1] = Vb m H N c ^ ( 3 .3) 

87T 

where the expression on the right-hand side is the four-dimensional result. 


4 Next-to-leading order 


4.1 Real emission contribution 

The real emission contribution to the differential decay width reads 

dr? = Kt/ - 


(4.1) 


dr B is divergent when the radiated gluon becomes unresolved (soft, or collinear with one of the 
b-quarks). In order to regularize it, we subtract an approximate decay rate, 


dI ^ Al = > (4-2) 

where the counterterm for processes with to + 1 partons in the final state is given by [32, 33], 

ra+1 

aix®iI 2 = E 

r= 1 



In eq. (4.3) the functions and 6>r°' 0) appearing in the right-hand side correspond to countert¬ 
erms which regularize the Pi\\p r collinear limit and the p r —> 0 soft limit. In order to avoid double 
counting in the overlapping soft-collinear region, we must add back a soft-collinear counterterm, 
C ir Sr°’° . The precise definitions of these subtractions are given in refs. [32, 33]. In our convention 
the indices of C^' 0 ' 1 are not ordered, C^’ 0 ' 1 = C^’°\ Since the sums over i and r in eq. (4.3) are 
likewise not ordered, the factor of | assures that we count each collinear limit precisely once. Fi¬ 
nally, the superscript (^ 1 ,^ 2 ) means that the corresponding counterterm involves the product (in 
colour or spin space) of an fd-loop unresolved kernel (an Altarelli Parisi splitting function or a soft 
eikonal current) with an ^ 2 -loop squared matrix element. Thus, (0, 0) means that we consider a tree 
level collinear or soft function acting on a tree level reduced matrix element. Such superscripts will 
appear also for other counterterms throughout the paper. For definitiveness, we spell out eq. (4.3) 
explicitly for H —► bb g (to = 2) below, 
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where the b, b and gluon carry the labels 1, 2 and 3. 

With the counterterms given in refs. [32, 33] it is straightforward to check that the difference 

dr NL° = dr R j 3 _ ^R.A, ( 4 .5) 


is integrable in all kinematic limits. Then, the regularized real contribution to the decay rate, 


r 3 LC V] = [ dr 3 L °] £=0 (4-6) 

is finite in four dimensions for any infrared-safe observable. An explicit calculation for the contri¬ 
bution to the total decay width from the real emission part plus subtractions yields 


rN LO [J= l] =T lu — C f 


,lo 1729 


450 


4.2 Virtual contribution 

The virtual contribution to the differential decay width reads 


dry = 


1 




2mjf'* T ‘ v "' H '""' bb ' b6/ ’ 
and is of course divergent in four dimensions. Its e-expansion reads (see eq. (A.2)) 

2tt <?ms \ m 2 H/ 


dr y = dr B ^ 7 5=f-^ ) C F 


-- - 2 + n 2 * +3L + 0(e) 


(4.7) 


(4.8) 


(4.9) 


where we have introduced the abbreviation L = lnf^j-J. In eq. (4.9), dr B denotes the d- 
dimensional Born decay rate as given in eq. (3.1). 

By the KLN theorem, the integral of the approximate decay rate precisely cancels the diver¬ 
gences of the virtual piece, so adding back what we have subtracted from the real correction, the 
virtual contribution becomes finite as well. We have performed the integration of the various sub¬ 
traction terms analytically in ref. [32] and here we only quote the result, which can be written 
as, 

J dr*+i = dr B ® if } (M m ; e), (4.10) 

where the ® product is defined in eq. (2.2) and the insertion operator is in general given by [32] 2 


J i 0) (Mm;e) 
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k^i 


(4.11) 


The variables yiQ and Yik,Q were defined in eqs. (2.8) and (2.9) and Q M is the total incoming momen¬ 
tum. The functions Cf \{yiQ\e) and sf {Yik,Q\ e) have been computed as Laurent expansions 
in e in ref. [32] and are recalled here up to finite terms in appendix B. We mention that there is no 
one-to-one correspondence between the unintegrated subtraction terms in eq. (4.3) and the kine¬ 
matic functions that appear in eq. (4.11). The latter are obtained from the former after summing 
over all unobserved quantum numbers (colour and flavour) in addition to integrating over the un¬ 
resolved momentum, and organizing the result in colour and flavour space. Loosely speaking, the 
integrated form of C,-f enters cf] and that of <sf^ enters sf However, we are free to assign 

2 The expansion parameter in ref. [32] was chosen implicitly, with the harmless factor l/fiMS suppressed. 

For the sake of clarity we reinstate the factor 1/Sy s here, as well as in all other insertion operators in eqns. (5.30), 

(5.34), (5.39) and (5.43) below. 
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the integrated form of C ir Sr°^ to either of the integrated counterterms and this final organization 
was performed differently in ref. [32] and in this paper. In ref. [32], the integrated form of C ir Sr 
was grouped into while here we find it more convenient to group it into C^ 0 ]. 

For H —> bb, with only two partons in the final state the colour connections factorize completely, 

T x T 2 = -C F . (4.12) 

Furthermore, momentum conservation implies that 


ViQ — U2Q — Y 12 ,q — 2/12 — 1 * 


(4.13) 


Thus, the insertion operator I 
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where, as indicated, we must evaluate all functions with arguments equal to one. The Laurent 
expansion of eq. (4.14) in e is, 
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where, for future reference, we have also provided the 0(e) part in terms of rational numbers and 
known transcendental constants. The uncertainty of the 0(e 2 ) numerical result, as well as those of 
all other numerical results we show affect the last quoted digit, unless specifically stated otherwise. 
It is easy to check that the expression 


dr 


NLO _ 
2 — 



dr 


R.Aj. 


J 2 


(4.16) 


is free of e-poles. Hence 

F^ L °[./] = ^ [dF^ LO ] e=0 (4.17) 

is finite in four dimensions for any infrared-safe observable. For the contribution to the total width 
from the virtual part plus integrated subtractions we find 

r^ LO [j = 1] = r LO ^ + \ c ^ • ( 4 - 18 ) 

Combining eqs. (4.7) and (4.18), we obtain the full NLO correction to the total decay rate, 

r NLO = r^ LO [j = 1] + r^ LO [j = 1] =r LO + \ c ^L S j ■ (4.19) 

As Cf = | in the conventions used, we recover the well-known NLO result [78-80]. 

5 Next-to-next-to-leading order 

5.1 Double real emission contribution 

The double real emission contribution to the differential decay width is 

dr 4 R = + E 1 M bij 2 + P)2l-^Sbbl 2 ) > ( 5 - 4 ) 
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and its integral over the phase space is divergent in four dimensions due to kinematic singularities 
emerging in unresolved regions. In order to regularize the singularities of eq. (5.1) due to two 
unresolved partons, we subtract an approximate decay rate, 


dr R R , A2 


2 mu 
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where the double unresolved counterterm for processes with m + 2 partons in the final state is [33] 
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(5.3) 


In eq. (5.3), the functions C^f \ c£+, GS+ 1 -* and 5^’ 0> denote counterterms which regularize the 
Pi I \Pr | IPs triple collinear, the Pi\\p r , Pj\\Ps double collinear, the pi\\p r , p s —> 0 one collinear, one 
soft (collinear+soft) and the p r —> 0, p s —> 0 double soft limits. The rest of the counterterms 
which appear in eq. (5.3) account for the double or triple overlap of limits, their role is to make 
sure that no multiple subtractions are performed in overlapping double unresolved regions. Thus, 
for instance, C irs CS^!^ accounts for the triple collinear limit of the collinear+soft counterterm, and 
the rest of the counterterms have a similar interpretation as suggested by the notation. The precise 
definitions of all functions appearing in eq. (5.3) were given in ref. [33]. As in our convention the 
collinear indices of counterterms and the sums over them in eq. (5.3) are not ordered, the factors of 
|, |, etc., are needed so that each limit is counted precisely once. 

After subtracting the double unresolved approximate cross section, the difference 

dr]+ - dr R R ’ A2 (5.4) 


is however still singular in the single unresolved regions of phase space. To regularize it, we also 
subtract 


dT 
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where A\ has been defined in eq. (4.3). To avoid double subtraction in overlapping single and double 
unresolved regions of phase space, we must also consider 


dr 
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The general formula for the iterated single unresolved counterterm is 
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where the three terms above are given by [33], 
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The interpretation of the various terms in eqs. (5.8)-(5.10) are suggested by the notation: for 
instance, C kt C^ t ’^ in eq. (5.8) accounts for the p k \\pt single collinear limit of the triple collinear 
counterterm, while, for example, S t C^f^ in eq. (5.9) represents the counterterm appropriate to 
the p t —¥ 0 soft limit of C^fK Thus, Ai 2 \M^ +2 \ 2 cancels the single unresolved singularities 
of the double unresolved subtraction term A 2 \M^ + 2 \ 2 ■ However, very importantly, it can also 
be shown [33] that A\ 2 \M.^ +2 \ 2 simultaneously cancels the double unresolved singularities of the 
single unresolved subtraction term Ai\M^ +2 \ 2 an d so properly accounts for the overlap of single 
and double unresolved subtractions. All of the counterterms appearing in eqs. (5.8)-(5.10) were 
precisely defined in ref. [33]. As before, the collinear indices and sums over them in eqs. (5.7)-(5.10) 
are not ordered, hence the appearance of the factors of l at various instances. 

With these definitions, the difference 

dr N N L° = dr RRj 4 _ dr f RA2 j 2 - drf R,Ai j 3 + dr RR,Al v 2 ( 5 . 11 ) 

can be shown to be integrable in all kinematic limits [33]. Thus, the regularized double real 
contribution to the decay rate 

r]+ LO [j] = J [dr^ NLO ] e=0 (5.12) 

is finite in four dimensions for any infrared-safe observable and can be computed with standard 
numerical techniques. For the total cross section (J = 1) at /i = win (L = 0) we find, 

r^ NLO [j = 1 ] = r LO 2 1.05(1). (5.13) 

This numerical value has been obtained by implementing eq. (5.12) in a fully differential parton level 
Monte Carlo program using four dimensional double real emission matrix elements and phase space. 
However, we have also reproduced the result by integrating the matrix elements and subtraction 
terms directly in d dimensions and then summing the separate contributions. We stress that this is 
a highly non-trivial cross check, as both calculations are very different conceptually and technically. 
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5.2 Real—virtual contribution 

The real-virtual contribution to the differential decay rate reads 


dlt v = 


1 




(5.14) 


which contains explicit e-poles coming from the one-loop matrix element and furthermore it is di¬ 
vergent in phase space regions where the gluon becomes unresolved. The explicit poles are cancelled 
by the integral of the single unresolved subtraction term in the double real emission contribution 
to the full NNLO decay rate, 


^ dT^’ Al = dr* 0 4 0) ( Pl , p 2 , p 3 ; e), 


(5.15) 


where the real emission differential decay rate is dT^ is given by eq. (4.1), while the insertion 
operator l[ 0> (p^,p 2 , pp e) is given by eq. (4.11). As there are only three partons in the final state, 
the colour connections that appear in the generic case in eq. (4.11) factorize completely, 


Thus, 

l[ 0) {PliP2,P3^) = 


T{T 2 = 


S r 


C A - 2C f 


P 

2-7T gMS \m 2 H 

C A 


C F 


and TiT 3 = T 2 T 2 = . 


C?J(yio; e) + C^(y 2Q ; e) - 2S[ 0) ^ 2 \y 12 , q] e) 


(5.16) 


C^(y 3Q -,e) + s^’ (1 ’ 2) (y L2j Q; e) - S^ 1 ’ 3 ^; e) - S^ 2 ’ 3 >(y 23>Q ;e) 

(5.17) 

Using the expressions in appendix B, it is straightforward to check that 




&s Se f P 

2tt S'MS \m 2 H 


2 Cp + C A 


+ "ttUa + 3Uf — — nfTn 
6 3 


(C A ~ 2Cf) In j/i 2 - CaO^Viz + lny 23 ) 
0(e°) 


hence the combination 


d T™ 


dr 


RR.Ai 


(5.18) 

(5.19) 


is finite in e. 

Nevertheless, eq. (5.19) is still singular in the single unresolved regions of phase space and 
requires regularization. We achieve this by subtracting two suitably defined approximate decay 

rates, dr 3 V ’ Al and ^/ 1 dr RR ’ Al ^ . First, we consider 


dT 


RV,A X 


2 mu 




(5.20) 


which matches the kinematic singularity structure of dr 3 v . The general definition of the real-virtual 
counterterm is [34], 


m+1 


r—1 

m+1 

+ E 


m+l 


m+1 


E + s< w > - E C .A° 


d) 


r =1 


. i =1 
i^r 

m+1 


E + 


(i,o) 


i= 1 
i^r 

m+1 


^’ 0) - Y, +E 


o) 


i= 1 
i^r 


i= 1 
iy^r 


(5.21) 
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The basic structure of this subtraction in terms of unresolved limits is the same as the tree level single 
unresolved counterterm in eq. (4.3). However, in accordance with the form of infrared factorization 
of one-loop QCD matrix elements [81-84], in eq. (5.21) we have terms with tree level collinear or 
soft functions multiplying (in colour or spin space) one-loop matrix elements (those with the (0,1) 
superscript), as well as terms with one-loop collinear or soft functions multiplying tree level matrix 
elements (denoted with the (1,0) superscript). The precise definitions of the functions appearing 
in eq. (5.21) are given in ref. [34]. 

Then we consider the counterterm, 

drr Al ) Al = ^d, (5.22) 

which matches the kinematic singularity structure of J] dr^ R,Al . In general, the counterterm is 
given by [34], 

m+1 

+iAC>i a ®'?”) = £ 

r =1 

m+1 

+ E 

r= 1 

i^r i^r 

The organization of this subtraction in terms of unresolved limits is again identical to the tree level 
single unresolved counterterm in eq. (4.3). However, for each limit, we have two types of terms, 
labeled by the different superscripts. The reason is as follows. This counterterm is built from 
the infrared factorization formulae for the product of a QCD squared matrix element times the 
I^ r> insertion operator of eq. (4.11). It turns out that these factorization formulae are sums of two 
pieces. Both of these involve the product of a tree level collinear or soft function times a tree level 
matrix element, but one piece is further multiplied by the insertion operator appropriate to 
the reduced matrix element, while the other is multiplied with a well-defined remainder function 
R [34]. Hence the superscripts on the various terms in eq. (5.23). 

It can be shown that the combination 




dP 


iNNLO __ 


dr 


RV 


dr 


RR,A X 


J 3 ~ 


dr 


RV.A! 


dr 


RR,A X 


L Jl J L 

is both free of e-poles and integrable in all kinematically singular limits [34], 
real-virtual contribution to the decay rate 


•h 


(5.24) 


Thus, the regularized 


rr LO W = ^[drr LO ] e=0 (5.25) 

is finite and can be computed numerically in four dimensions for any infrared-safe observable. For 
the total cross section ( J = 1) at [i = m_ h ( L = 0) we find, 

r^ NLO [J = 1 ] = r LO ) 2 69.35(1). (5.26) 

As for the double real emission contribution, the numerical result of the Monte Carlo program in 
eq. (5.26) has been reproduced by integrating the real-virtual matrix element and the subtraction 
terms separately in d dimensions and summing the contributions. 


- 11 - 











5.3 Double virtual contribution 

The double virtual contribution to the differential decay rate reads 


dry v = 


2 m H 


d <j> 2 (m 


(5.27) 


which contains explicit e-poles coming from the two-loop matrix element and the square of the 
one-loop matrix element: 


dry v = dr B 


2Cf 


+ 


Og 

2tt S m 


8 7r 2 11 , 

9 + 12 + Y L]CaCf + 


11C a C f 


+ (6 + 4 L)Cl - ntT R C F 


17 


- 2tt 2 + 6 L + 4 L 2 Cl - 


+ ( ^-2^ 2 -14C 3 +4(2-^ 2 )L + 3L 2 + ^L a )C^ 


+11 + j L )”> T ^ 


- + 0(e") }. 


- + -L IrifT^CF 


(5.28) 


The finite part of dry v is also known exactly [85] which we recall in appendix A (see eqs. (A. 3) 
and (A. 4)). In order to regulate these poles, we add the integrals of the counterterms which have 
been subtracted in sections 5.1 and 5.2. The KLN theorem then ensures that, provided the physical 
observable we are to compute is infrared-safe and our subtraction scheme is internally consistent, 
the ensuing result will be free of infrared divergences. It is our task in this section to verify that 
this is indeed the case. 

Let us begin with the integral of the double unresolved subtraction term, eq. (5.2), which can 
be written as, 

f 2 dI t+ 2 A2 = dr B ® 4 0) ({p} m ; e), (5.29) 

where the insertion operator has five contributions according to the possible colour structures, 


4 0) (M;e) = 


a s S e ( /i : 


2tt sms y Q 2 


E 


C { °> (y iQ ; e) T 2 + £ C ( 2 % ( y lQ , y jQ , Y^q ; e) T 


3= 1 


E 


SF ,U, ) (Yji,Q-,e) C a +Y, ^ 2 ,i {yiQ>Yij,Q>Y iljQ ,Yj liQ -, e) T 


3,1 =1 L 


TjT t 


m m 


E E 


i,k= 1, j,l= 1, 
k^i 


(Yik,Q,Y i j i Q,Y i i t Q,Yj k 'Q,Y k i i Q,Yji t Q- : e){T l T k ,T j Ti}j . 

(5.30) 


The kinematic functions in eq. (5.30) have been defined and computed as expansions in e in 
refs. [40, 41]. Again, there is no one-to-one correspondence between the unintegrated double un¬ 
resolved subtraction terms in eq. (5.3) and the kinematic functions that appear in eq. (5.30). The 
latter are obtained from the former after integration over unresolved momenta and summation over 
unobserved colours and flavours. This remark applies to the rest of the insertion operators to be 
discussed below. 
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For H —» bb, the colour connections that appear in eq. (5.30) are simply given by eq. (4.12), 
and the kinematic variables simplify as in eq. (4.13). Furthermore, when evaluating eq. (5.30) the 
coincidence of certain summation indices is allowed. In particular, i in the second line need not 
be distinct from j and l , while in the last line we only require that i and k as well as j and l 
are different, with no further restrictions, as shown in the formula. As a result, some indices of 
kinematic functions coincide once we explicitly write out eq. (5.30). Specifically, since in our case 
there are only two hard partons in the final state, only CS 2 °-’^’^ and g( 0 b(*,fc).(*.O a pp earj w hil e the 
more general functions or SOld,O,(TO are a t) Sen t from the sum, as those require at least 

three hard partons if all indices are different. In such cases we also simplify the list of arguments 
of the functions so that we do not display arguments that are the same or identically zero. For 
instance, in CS 2 °]’^’^ if * = j , then Yi^Q — 0 and Yu,q = Y)/,q. Hence, CS^ 0 ]’^’^ is a function of 
UiQ and Yh q only. Similarly g(°b(c 0 , 0,0 ^gpg^g j us t 0 n the variable Y^ q. Then, we obtain the 
^ 2 °^(PiiP 2 ; e) operator, 


-f2 0) (Pi>P2;e) = 


27r gMS \m 2 H 


2 Cl 


C 2 °n(l; <0 + c'l(l, 1,1; e) - 2CS®’ (1 ’ 2) (1,1; e) 


+ 4S (0),(1,2)(1,2) 


(l;e) 


-2C F C , AS^ ),(1 ’ 2) (l;e) 


(5.31) 


whose e-expansion is 


-f2 0) (Pi,P2;e) = 




27T S'MS \mjj 




2C f) ? 


/29CaCf „ 2 nfT R 6V \ 1 

V^2~ +6Cp 3~ J? 


68 _ 7tt2\ f 170 _ 8^\ 2 _ 14n f T R C F 

9 12 J A F + l 9 3 / F 9 


1 


301 

216 


37n2 r , f 6U9 

^- + t )CaC f +I — 


47tt 2 
~18~" 


- 70C 3 C$ 


+ I -g + ^ ] n,T R C r 


- 227.559 C A C F - 236.532+ 30.9273n f T R C' F + O(e) 


(5.32) 


The coefficients of the poles are all given in terms of rational numbers and known transcendental 
constants. 

Next, we consider the integral of the iterated single unresolved subtraction term, eq. (5.6), which 
can be written as, 


dr RR,A 12 =dr B g) J (0 ) ( Mm;e)) 


(5.33) 


where the insertion operator in general has the same structure in colour and flavour space as I 


(o) 

2 
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in eq. (5.30), 


J i2(W; e ) = 


a s S e ( p 


2ir S MS \Q 2 


E C S;<(»«?; £ )+ E c i 0 ifc(^> £ ) T 


*=i l 


fc=i 

k^i 


+ E s i°2 (J e )^A + E GsSi W^,q> ^,q; c)t: 


j,i=i 




m m 


E E si 0 2 ) ’ (l ’ fe)a ’ i) (y ifc ,Q,F u -Q,r Ii , o ,F Jfc ,Q,y fciiO ,F j; , Q ; e ){r i T fe ,r J T i } }>. 


i,/e=l j,l =1 


(5.34) 


The kinematic functions in eq. (5.34) have been defined and computed as expansions in e in ref. [39]. 
The discussion below eq. (5.30) applies to eq. (5.34) as well, hence, using eqs. (4.12) and (4.13), we 
obtain the (pi,P 2 ! e) operator, 


J i2 } (Pii^; e ) = 


Q s S e ( p 2 

2tt5'MS 


2 Cl 


Cg ? (l; e) + Cg„(l, 1,1; e) - 2CSgf’ 2) (l, 1; e) 
+ 4S (0),(l,2)(l,2) (1;e) l ^ 2CFC . A s(° 2 )-(b 2 ) ( l; e )' 


(5.35) 


whose e-expansion is 
I { i2(Pi,P2',e) = 


a s S e / /i 2 
2 tt V TO 1/ 
155 


4C 2 


CaCf , 10/02 2nfT R CF\ 1 

-^- + 12 G f 3 ^ 


18 


-T^r+r' C A C F + —- 


788 16tG 


25 


r 2 - 

°F 


31n f T R C F 


5911 101 In 2 49 tt 2 _ . _ _ 

^r + — + ^r + 4 2c 3 ]c A c P - 


/116497 296tt 2 

^ 4500 + 45 


'71 202 In 2 8 tt 2 \ 

+ l 36 9— + “9“ J ^fTRC'F 


+ 215.508CaC f - 717.881C| + 22.1494n f T R GV + O(e) \. 


104C 3 C 2 


(5.36) 


As in the case of I^\ the coefficients of the poles are all given in terms of rational numbers and 
known transcendental constants. 

Turning to the integral of the real-virtual single unresolved subtraction term, eq. (5.20), we 
find [35] 

[ dr^+f 1 = dr^ <g> 4 0) ({p} m ; e) + dT« ® e), (5.37) 

where the insertion operator 7]° is given in eq. (4.11), expanded to sufficiently high order in 
eq. (4.15) to obtain the first term on the right-hand side in eq. (5.37) to O(e), while the operator 
in general reads 


I { i\{p} m ;e) = l[^ B ({p} m -e) M m ;e). 


(5.38) 
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/ -j \ Q 

The unrenormalized operator I\ ’ has the following structure in colour and flavour space, 




a s S f ( n 


2tt S ms \Q 2 


E c [f(w £ ) c A r? + E s 


(1 


f Yik,Q\e) C A T,T k 


k= 1 
k^i 


m m 

EE S 

k =1 Z=1 
k^i l^i,k 


(1 


(Y ik>Q ,Y atQ , Y kliQ] e) E fabcT^T^ 


(5.39) 


The bare kinematic functions in eq. (5.39) have been defined and computed as expansions in e in 

/ l N g 

ref. [35]. Using eqs. (4.12) and (4.13), the unrenormalized I\ ’ operator becomes 


1 1 (Pi,P2;e) = 


a s / M 2 

27T5-MS 


2C a C f 




(5.40) 


The term involving triple colour correlations on the second line of eq. (5.39) does not contribute, 
the triple sum over i, k and l being empty because we cannot form a triplet of distinct indices. The 


e-expansion of the bare insertion operator l\ ’ reads 

2 


1 1 {Pi,P2;e) = 


a s S e / E_ 

2tt ^ms \m 2 H 


(i ),b 

l 

C A C F 3 6' a C'f 


2e 4 


2e 3 


H ^ VaCf + ( - ? + ^ )C? 


450 2 

(580571 43tt 2 


V 6750 30 

292.930C A C F + 134.720Cf + O(e) 


-15Ca)c A C F+ f-- ^ 


^2 


661 13184 In 2 71 tt 2 


45 


38^3 JC'p 


(5.41) 


Also here we see that the pole coefficients are all given in terms of rational numbers and known 
transcendental constants. 

Finally the iterated integral of the double real single unresolved subtraction term, eq. (5.22), 
can be written as, 


dT; 


RR.AO Ai 


1 x Jl 


m+2 


= dTf 


\ {4 0) (W m ; e), I?\{pU 0} + e) 


(5.42) 


where the insertion operator is given in eq. (4.11), expanded to sufficiently high order in 
eq. (4.15) to obtain the first term on the right-hand side of eq. (5.42) to O(e), while l\ 0 \ (}> in general 
reads 




a s S e (ji 


27T s MS \Q 2 


E C 1°M (y<o; e) Cat 2 + E si 0 i 0) (Y ik , Q - e ) c A t,t. 


k —1 
k^i 


(5.43) 

The kinematic functions in eq. (5.43) have been defined and computed as expansions in e in ref. [35]. 
Using eqs. (4.12) and (4.13), we obtain 


I i ) f\pi,P2-,e) = 


a s S e / 6 

2tt S'MS \m 2 H 


2 CaC f 


C^(l;e)-S)E’^(1U) 


( 0 , 0 ), ( 1 , 2 )/. 


(5.44) 
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whose e-expansion is 


A °i ) (Pi,P2-,e) = 


S f 


2n S^s \m 2 H 




e-i 2 


f C\Cf 2nfT R CF\ 1 
” I” 3~ 3 ) e 3 


-if + - 2 > c * c - 


5 - 


47r 2 


C£- 


31n f T R C F 


622583 101 In 2 502 tt 2 


3375 


/393797 13184 In 2 

+ + 


45 

274t r 2 


50C 3 C A C F 


+ 


v 13500 
/11557 


225 


45 


- 66 C 3 cl 


202 In 2 87 t 2 . 

—+ -g- ) n { T R C F 


V 2700 

- 15.2343CACF - 318.099C| + 46.4407n f T R C F + 0(e) 


(5.45) 


All the pole coefficients are again given in terms of rational numbers and known transcendental 
constants. 

Using eqs. (5.28), (5.32), (5.36), (5.41) and (5.45), it is straightforward to check that the 
regularized double virtual contribution 


dr NNLO = J dr yv 


ipRRjAo j-pRR,Ai 

dr 4 2 -dr 4 


'1 L 


dr 


RV.A! 


dr 


RR,Ai 


J 2 (5.46) 


is free of e-poles. Hence, the regularized double virtual contribution to the decay rate 


r NNLO [J]= ^[ dr NNLO 


(5.47) 


is finite for any infrared-safe observable and can be computed numerically in four dimensions. For 
the total cross section (J = 1) at /i = run (L = 0) we find, 


r^ NLO [j = 1 ] = -r LO (yy) 2 41.25(1). 


(5.48) 


We note that the error estimate of the above result comes entirely from the uncertainty associ¬ 
ated with the numerical computation of the finite parts of the insertion operators. The statistical 
uncertainty of the Monte Carlo integration over the two-parton phase space is completely negligible. 
Finally, summing eqs. (5.13), (5.26) and (5.48), we obtain 

pNNLOjjr _ — pNNLOjqj _j_ pNNLOjqj _j_ pNNLO |q] _ pLO ^^29.15(2), (5.49) 


to be compared with the know analytic result 

F NNLO [J= 1] =F LO 29.146714... . 


(5.50) 


6 Inclusive and differential results 

In this section, we show that using the fully differential two-, three- and four-parton contributions 
of eqs. (3.1), (4.5), (4.16), (5.11), (5.24) and (5.46), we can make predictions for any infrared-safe 
jet cross section with jet functions J n (n = 2,3 and 4) defined in d = 4 dimensions. 
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/V m H 


Figure 1 . Scale dependence of the inclusive decay rate at LO, NLO and NNLO accuracy. The estimated 
uncertainty on the numerical results is too small to be appreciated. 


The inclusive decay rate is obtained by setting J = 1 and is given by the sum of the leading 
order width (3.3) and the NLO (4.19) and NNLO (5.49) corrections. At = m# we obtain 


p pLO 

1 NNLO — J- 


Os 17 
7r 3 


(^) 2 29.15(2) , 


( 6 . 1 ) 


in agreement with the known analytic prediction [78-80]. In figure 1, we compute the inclusive 
decay rate at /r = m#/2 and /i = 2m h and compare it to the known analytic result for the scale 
dependence, finding excellent agreement. 

To illustrate the impact of NNLO QCD corrections on differential distributions, we apply the 
Durham jet algorithm [86] with resolution parameter y cut = 0.05 to cluster final state partons and 
order the resulting jets in energy. In the top panel of figure 2 we show the energy distribution of 
the leading jet in the rest frame of the decaying Higgs boson for two-jet events. In ref. [47] the same 
distribution was computed for jets clustered according to the JADE algorithm with y cut = 0.1. We 
have repeated that calculation and found excellent agreement with the published results. However, 
for two-parton kinematics the energy of the leading jet is just E max = mn/2., so at leading order 
the leading jet energy distribution is a delta function. Furthermore, double unresolved subtractions 
for four parton matrix elements, as well as single unresolved subtractions for three parton matrix 
elements also contribute to this distribution only at i? max = mu /2. Then, to show the subtraction 
method at work on an observable that has a non-trivial distribution already at leading order, we 
consider the absolute value of the pseudorapidity of the leading jet, \rji |, with respect to an arbitrary 
axis. The effect of higher order corrections on this distribution is shown on the bottom panel of 
figure 2. In this last illustrative example we note that going from the leading order to NNLO, the 
uncertainty bands shrink, and that the NNLO band falls within the NLO band, thereby showing 
the good convergence of the perturbative series. 
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Durham clustering at y cut = 0.05 



^max / r/ > / / 



0.0 0.5 1.0 1.5 2.0 2.5 3.0 


M 


Figure 2. The plots show the normalized distribution of the leading jet energy i? max (top) and the 
distribution of the absolute value of the pseudorapidity 1 771 1 of the highest energy jet (bottom) at LO, NLO 
and NNLO accuracy. The bands show the dependence on the renormalization scale corresponding to the 
range y £ [mjj/2, 2mjj]. Jets have been clustered using the Durham algorithm, the resolution parameter 
for jet clustering was set to y cu t = 0.05. 


The bands in both distributions in figure 2 correspond to the envelope of varying the renormal¬ 
ization scale in the range £ [ran/ 2 , 2 mjy]. 
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7 Conclusions 


In this paper, we have computed the fully differential decay rate of the SM Higgs boson into b- 
quarks at NNLO accuracy in a s , by implementing a general subtraction scheme developed in a 
series of papers for the computation of QCD jet cross sections at NNLO accuracy [31-41]. 

We have shown that our subtractions render both the double real and real-virtual contributions 
to the NNLO correction integrable in four dimensions. We have also presented the integrated forms 
of our subtraction terms with pole coefficients evaluated analytically, while the finite parts were 
given numerically. We confirmed that the sum of the double virtual contribution and the integrated 
subtractions is free of infrared singularities as required by the KLN theorem. We have implemented 
our computation in a parton level Monte Carlo program and presented illustrative examples of 
differential distributions at NNLO. 

The successful application of our subtraction scheme reported here opens the way to the com¬ 
putation of other, more involved processes and is also encouraging to further developments of the 
scheme to deal with initial state radiation. These directions of development are under way and will 
be the subject of further publications. 
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A Matrix elements 


We present the matrix elements in the form used in our parton level Monte Carlo program. In 
particular, in our scheme we need the the four-parton tree level and the three-parton one-loop 
matrix elements only up to finite terms in e. Higher order terms must of course be included when 
integrating the matrix elements and subtraction terms separately in d dimensions. When needed 
for our cross checks, we take these higher order terms directly from ref. [47]. 


A.l Two partons 

For H —> bb at tree level we have 


l-^bsl 2 = 2 vl m 2 H N c 


(A.l) 


We computed the one-loop correction and obtained 


2SK?|A-®=“- * 


2tt S ms V rn j-i 


KIWf! - 4---2- 
3 


3 L 


1 


-( 4 +^- 4 Cs + ^ 2 )e-( 8 - 7 r 2 + C3 + ^-^-^ 3 )e 2 


0(e 


(A.2) 


We used the formula at two loops as given in ref. [47] : 




Os 


A 


27r SMS y mjj 
7 r 2 

+ 12 


2 /IICaCf 
4 + y 4 


4“ 3Cf — rifT R C F ) — 


- ^L) C a C f + I ^ - 2tt 2 - 3L ) C 2 - ( ^ - ^L) 


nfT R CF 


961 13Cs 
53 3 tt 2 


S3 _ 11 L 11 L 2 j c c 
216 2 2 6 / 


4 -^L+-L 2 \C F + 


^+2L-^L 2 )mT R C F 


467 733tt 2 

162 + 216 


1 7 _ 557T 2 _ 20 ^ 3 + 


92Cs 11 tt 4 


24 


9 360 

43t r 4 / 9 
90 ~~ l 4 


53 

12 

57T 2 

AT 


5571"* 

”36" 


l + y l2 ~^ l3 ) CaCf 




^-^-^-fi + ^L-2L 2 + ^)n f T R C F 

81 54 9 V 3 9 / 9 ' 


0(e) 


(A.3) 


We checked that the poles of this expression satisfy the general formula given in ref. [87], while the 
finite part agrees with that in ref. [85]. The square of the one-loop matrix element is 


25R(A4«|A4«) = 


S f 


27r S MS \ TO 7f 






, 5t r 2 _ 9 3 2 \ 1 
7- — -4<; 3 --L+-L )- 
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A.2 Three partons 

For H —> bb g at tree level we have 

^S/ = 8 -iBs" 2 ‘^Si 2 cp4 
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(1 — e )y 23 , (1 -e)y 


13 


y 13 


2/23 


2 yi 2 
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+ 2 — 2e 


(A.5) 


At one loop, we use the e-expansion of the formula from ref. [47], which we checked numerically 
against GoSam [88 , 89], 
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where 


i?(x, 2 /) = Li 2 (1 - x) + Li 2 (1 - y) + lnxlny - — . 
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(A.7) 


A.3 Four partons 

In our computation we need the H —> four partons squared matrix elements at tree level in d, = 4 
dimensions. We checked our formulae, presented below, with GoSam [88 , 89]. 

For H —> bbgg we have 


i-^SL-i 2 = (s™^ 2e ) 
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where 
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where 
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where 


A bbgg(Pl^P2,P3,P4.) 
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B J, 0 ' 1 insertion operator to 0(e) 

We present the /)°^({p} m ;e) insertion operator in eq. (4.11) to 0(e). More precisely, we give the 
e-expansion of the kinematic functions C^ 0 - (x, e) and S w ] =t i c h. appear in eq. (4.11) up to 
and including finite terms. 

Starting with C^ 0 ,- (x, e), we have 

C^(x, e) = [Ci 0 ) ] ? 9 (x, e) - [C ir S<°>](e) , (B.l) 

cQ(x,e) = ^Ci 0 ) ] gs (x,e) +n f [Ci 0 ) ] gg -(x,e) - [C fr S<°>](e) , (B.2) 
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